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Inspired by the recent developments in modeling and analysis of reaction networks, we provide a geometric formulation of the 
reversible reaction networks under the influence of diffusion. Using the graph knowledge of the underlying reaction network, 
the obtained reaction-diffusion system is a distributed-parameter port-Hamiltonian system on a compact spatial domain. 
Motivated by the need for computer based design, we offer a spatially consistent discretization of the PDE system and, in a 
systematic manner, recover a compartmental ODE model on a simplicial triangulation of the spatial domain. Exploring the 
properties of a balanced weighted Laplacian matrix of the reaction network and the Laplacian of the simplicial complex, we 
characterize the space of equilibrium points and provide a simple stability analysis on the state space modulo the space of 
equilibrium points. The paper rules out the possibility of the persistence of spatial patterns for the compartmental balanced 
reaction-diffusion networks. 
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1 Introduction 

Reaction-diffusion systems model the evolution of the 
constituents distributed in spac e under the inf l uence of 



chemical reactions and diffusion Smoller (1994); Temam 

f 97| ). These spatially distributed models are essential 
The understanding of many important phenomena 
concerning the development of organis ms, coordinate d 
cell behavior, and pattern formation Murray (2003). 



Guided by the models of reaction-diffusion equations, 
designing multicellular systems for pattern formation is 
one of the present research topics in synthetic biology, 
with application foreseen in tiss ue engineer i ng, bi oma- 
terial fabrication and biosensing Basu et al (2005). 



The mathematical model of a chemical reaction that 
proceeds in an ideally diluted environment is a dynam- 
ical system x = f(x). A physicochemical interpretation 
would be that of a reaction system involving m species 
xi, . . . under the influence of kinetics modeled by 
the vector field /. If a well- mixed hypothesis is not rea- 
sonable, a more appropriate model is that of reaction- 
diffusion equations 

Ox 

— = div(D(x)gradx) + f(x) , (1) 
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withx:= (£i(£,t),...,£ m (£,f)) T : (M,R+) -> M m , and 
a positive semi-definite diagonal diffusion matrix. The 
operators grad and div act component-wise with respect 
to the local coordinates £ = (£i, . . . , £ n ) of the compact 
spatial domain M C W 1 . The constraints acting on the 
system from the outside M impose appropriate bound- 
ary conditions, which together with all the other techni- 
calities will be discussed later. 

Motivated by the rece nt advances in the network c on- 
trol and graph theory, van der Schaft et al (2012) of- 
fers an elegant for mulation for the dynamic s of reversible 



chemical reactions Horn fc Jackson] (|1972[); Horn (1972) 



and 

the chemical reaction networks considered in 



Feinberg| ( |1987[ 1995). The graph descri ption of 



derl 



Schaft et al\ ( |2012[ ) has a direct thermo dynamical inter- 
pretation, which can be regarded as a graph-theoretic 
version of the formulation derived in the work of [Oster 
fc Perelson|fll974D;|Oster et al .|(|1973D. Based on this for- 
mulation, van der Schaft et al. |2012 ) characterizes the 
space of equilibrium points and provides a dynamical 
analysis on the state space modulo the space of equilib- 
rium points. 

After a brief review of balanced reaction networks closely 
following the exposition of |van der Schaft et al. ( |2012[ ), 
we introduce a Dirac structure that captures the geome- 
try of reaction-diffusion systems. We start from the fact 
that the considered reaction systems are defined with 
respect to a finite Dirac structure on a manifold. This 
means that the reaction system from a network modeling 
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perspective can be described by a set of energy-storing 
elements, a set of energy-dissipating (resistive) elements, 
and a set of ports (by which the interconnection is mod- 
eled) , all interconnected by a p ower-conserving intercon- 
nections | van der Schaft| ( 2000 ) . 



From a control and interconnection viewpoint a prime 
desideratum is to formulate reaction diffusion systems 
with varying boundary conditions in order to allow 
energy flow through the boundary, since the interac- 
tion with the environment takes the place through the 
boundary. The Stokes-Dirac structure offers a geomet- 
ric framework for this and allows us to model balanced 
reaction-diffusion systems as distributed-parameter 
port-Hamiltonian systems. 

It is well-known that adding diffusion to the reaction sys- 
tem can generate behaviors absent in the ode case. This 
primarily pertains to the problem of diffusion-driven in- 
stability which constitutes the basis of Turing's mecha- 
nism for patte rn formation Turing ( 1952 ), Nicolis & Pri- 



gogme 



(1977). Here, the port-Hamiltonian perspective 



2 Balanced Chemical Reaction Networks 

Stoichiometry. Consider a chemical reaction network 
involving m chemical species (metabolites), among 
which r chemical reactions take place. The basic struc- 
ture underlying the dynamics of the vector x G R+ of 
concentrations o^, i = 1, . . . , m, of the chemical species 
is given by the balance laws: x = Sv, where S is an 
m x r matrix, called the stoichiometric matrix. The 
elements of the vector v G R r are commonly called the 
(reaction) fluxes. The stoichiometric matrix 5, which 
consists of (positive and negative) integer elements, 
captures the basic conservation laws of the reactions. 

The Complex Graph. The network structure of a 
chemical reaction network cannot be directly captured 
by a graph involving the chemical species, because, in 
general, there are more than two species involved in a 
re action. Following the app r oach o r iginat ing in the work 
of |Horn fc Jacksonj ( [19721 ); |Horn| fll972| ) and |Feinberg 
3l) , we will introduce the space of complexes. 



1995) 



( [19871 

The set of complexes of a chemical reaction network is 



permits us to draw immediately some conclusions re- 
garding passivity of reaction-diffusion systems, but also 
to claim the spatial uniformity of the asymptotic behav- 
ior of balanced reaction-diffusion system. 

In the second part of the paper, by adopting a discrete 
differential geometry-based approach and discretizing 
the reaction-diffusion system in port-Hamiltonian form, 
apart from preserving a geometric structure, a compart- 
mental model analogous to the standa rd one is obtained 



Jacquez ( 1972 ); Jovanovic et al. ( 2008 ). Furthermore, we 
show the asymptotic stability of compartment al model 
and verify this result on an example of glycolisis path- 
way reactions. 

Notation. The space of n dimensional real vectors con- 
sisting of all strictly positive entries is denoted by R™ 
and the space of n dimensional real vectors consisting 
of all nonnegative entries is denoted by WL. l m denotes 



simply defined as the union of all the different left- and 
right-hand sides (substrates and products) of the reac- 
tions in the network. The expression of the complexes in 
terms of the chemical species is formalized by an m x c 
matrix Z, whose p-th column captures the expression of 
the p-th complex in the m chemical species. 

Since the complexes are left- and right-hand sides of the 
reactions they can be naturally associated with the ver- 
tices of a directed graph, with edges corresponding to 
the reactions. The resulting graph is called the complex 
graph. The complex graph is defined by its c x r inci- 
dence matrix B, c being the number of vertices and 
r being the number of edges. There is a close relation 
between the matrix Z and the stoichiometric matrix S, 
which is expressed as S = ZB. For this reason we will 
call Z the complex stoichiometric matrix. Hence, the 
relation x = Sv can be also written as x = ZBv, with 
the vector Bv belonging to the space of complexes R c . 



a vector of dimension m with all entries equal to 1. Balanced Mass Action Kinetics. A vector of con- 



The time-derivative ^(t) of a vector x depending on 
time t will be usually denoted by x. Define the map- 



ping Ln : 



— » 



x 1 y Ln(x), as the mapping 



whose z-th component is given as (Ln(x)) i := ln(^). 
Similarly, define the mapping Exp : R m — >• R+, x \-> 
Exp(x), as the mapping whose z-th component is given 
as (Exp(x))^ := exp(x^). Also, define for any vectors 
x,z G R m the vector x • z G R m as the element- wise 
product (x • z) i := XiZi, i = 1, 2, . . . , m, and the vector 
| G R m as the element- wise quotient (| ) := i = 
1, • • • , m. Note that with these notations Exp (a; + z) = 
Exp (x) -Exp (z) andLn(x-z) = Ln(x)+Ln(z), Ln (|) = 
Ln(x) — Ln(z). On an n-dimensional smooth manifold 
M, the space of smooth scalar valued functions will be 
denoted by C°°(M), while the space of vector valued 
functions will be C°°(M;R n ). We employ the notion 
C™(M) := C°°(M) x • • • x C°°(M), with the product 
being taken m times. 



centrations x* G R+ is called an equilibrium for the dy- 
namics x = Sv(x) if Sv(x*) = 0. Furthermore, x* G R+ 
is called a thermodynamic equilibrium if v(x*) = 0. 
Clearly, any thermodynamic equilibrium is an equilib- 
rium, but not necessarily the other way around (since in 
general S = ZB is not injective). 

Consider the j-th reaction from substrate Sj to product 
Vj , described by the mass action rate equation 

Vj{x) = fcf w exp {Z}.Ln(x)) - fc; ev exp (Z^.Ln(x)), 

where Z$ 6 and Z*p. denote the columns of the complex 
stoichiometry matrix Z corresponding to the substrate 
and the product complexes of the j-th reaction. Then 
x* G Rip is a thermodynamic equilibrium, i.e., v(x*) = 

0, if and only if Kj(x*) := kf rw exp (z^Ln(a?*)) = 
kj ev exp (^Zj>.Ln(x*)^ , for j = 1, • • • , r. The mass action 



2 



reaction rate of the j-th reaction now can be written as 

v(x) = -/C(z*)£? T Exp (z T Ln (^)) , 

where JC(x*) is the r x r is a positive diagonal ma- 
trix of balanced reaction constants given as JC(x*) := 
diag(^iO*),--- ,K, r (x*)). 

The dynamics of a balanced reaction network takes 
the form 

x = -ZB)C(x*)B T Exp (Vim . (2) 

This form will be the starting point for the analysis of 
balanced chemical reaction networks in the rest of this 
paper. Furthermore, we shall assume the validity of the 
global persistency conjecture, which states that for a pos- 
itive initial condition xo G R+, the solution x of (2) sat- 
isfies: liminf t ^ 00 x(t) > 0. The global persistency con- 
jecture re^entlvwasjroven for the single linkage class 



Anderson] ( |2Q 1 1 1 ) , but for the system (2) remains 



an open problem. 

Stability of Balanced Reaction Networks. It follows 
that once a thermodynamic equilibrium is given, the 
set of all thermodynamic equilibria is described by the 
following proposition. 



a thermodynamic equilibrium, then the set 
of all thermodynamic equilibria is given by 



£ := {a 



S T Ln(x**) = S T Ln(x*)}. (3) 



Making use of the formulation of the dynamics of bal 



anced reaction networks in (2), in van der Schaft et al. 



(2012) it was shown that all equilibria of a balanced re- 



action network are actually thermodynamic equilibria, 
and thus given by (3). A similar result was obtained in 
the classical papers Horn (1972 );|Hornfc Jackson] ( |1972[ ); 
Feinbergj ( |1995| ) , but also m ]Sontag| ( 12001 p , for a dlrTer- 
ent class of chemical reaction networks! 



Theorem 2 ( |van der Schaft et aE1 ( |2012D ) Consider 
a balanced mass action reaction network given by (2), 
and in addition assume the global persistency property, 
then for every initial condition x(0) G M.™, the species 
concentration x converges for t — ^ oo to E. 

In the remaining of the paper, we will extend some of 
the results of balanced chemical reaction networks to 
spatially distributed systems. 

3 Geometric Formulation 

In classical field theories the geometric content of the 
physical variables is usually expressed by identif ying 
them with differential forms of appropriate order |van| 




ifb,e b ) 



fay) 

Fig. 1. Reaction-diffusion system as a dissipative distributed 
port-Hamiltonian system. The conjugate variables u and y 
represent the inflows and the outflows of the reaction dy- 
namics. In this paper the reaction system is considered to 
be closed; that is, either u = or y = 0. 



Proposition 1 (van der Schaft et al. ( 2012| )) Let boundarv 



der Scha ft fc Maschke| (|2002). Previously in |Seslija et 
^ Q2010P we have offered a Stokes-Dirac structure 
lor reaction-diffusion systems. In this paper we avoid 
the exterior formulation and introduce a Dirac struc- 
ture for reaction-diffusion systems defined on the space 
of smooth functions on a Riemannian manifold with 



The formulation of a reaction-diffusion system as a port- 
Hamiltonian system on a compact n-dimensional smooth 
Riemannian manifold M with boundary dM is given 
as follows. We identify the mass density variables with 
an m component vector of scalar valued functions, that 
is x G Cm(M). The influence of the external world 
(reaction-diffusion system) to the system (outside world) 
is modeled through the boundary efforts e^ G (dM) 
(and boundary flows fa G C^(dM)). The reaction part 
is in its nature finite-dimensional and as such is mod- 
eled as the interconnection of the atomic elements, each 
of them characterized by a particular energetic behavior 
(energy storing, energy conversion or dissipation) . Each 
of these elements can interact with the environment by 
means of a port — a couple of inputs and outputs whose 
combination gives the power flow. The transport of the 
constituents in space is governed by the laws of diffusion, 
which is modeled as a thermal damping by termination 
of the appropriate ports (see Figure 1). 

Let the space of flows be T and its dual the space of 
efforts be £ . We set T := T x ® J~d © T r ® T\> and 
£ := £ x © E d © E r © £ b , where T x = E x = C%(M) is the 
carrier space of concentrations, Td = Ed = C^(M;R n ) 
the space of gradients, T r = E r = C^(M) the reac- 
tion flows, Th = C^(dM) being the space of boundary 
fluxes, and Eh = C^(dM) be the space of chemical po- 
tentials restricted to the boundary dM. 



3 



The inner product in C^(M) is given by 



4 Reaction-Diffusion Dynamics 



A I 



where d£ := d£i • • • d£ n is the volume element on M. 
Similarly, the inner product on the boundary dM is de- 
fined by 



JdM 



where d^A : = d£i * * * d£ n _i is the volume form on the 
boundary dM. Analogously, we defined the inner prod- 
uct in C%(M;R n ) and C%{dM\ W 1 ' 1 ). 
A non-degenerate pairing between T and £ is defined by 
the following bilinear form on T x £ with values in R 

i/ffl fl fl fl 1 1 1 1\ (f2 f2 f2 f 2 2 2 2 2\\\ 
\\\J x i J d ' J 7' J b 5 G x 5 °d? °r ' °6 / ' \Jx'>Jd'>Jr'>Jb ' °a; ' °d ' °r ' °fe / // 

:= (/a; J 6 x)lI i (M) *di e d) L^(M) ^ r ^ e r ) L^M) 

+ {fxi e x) L^(M) (fdi e d) L^(M) ^ r ^ e r ) L^M) 

+ (fbi e b) L 2 (dM) + {fbi e b)i 



f L* m {dM) 



where {f x , p d , ft, ft) e T and (4>4>4>4) G <M = 
1,2. 

The Stokes-Dirac structure that underpins the geome- 
try of reaction-diffusion systems is a maximally isotropic 
subspace of J 7 x £ . The following theorem gives the con- 
struction of a such Dirac structure. 

Theorem 3 Define V (Z T x £ by 

V : = {(fx,fdjr,fb,e x ,ed,e r ,eb) e J 7 x £ \ 



( fx ) 




/ div 


fd 




grad 


\e r ) 




\ Z T 




1- 


f tr 






I i/-tr 



(4) 



U/ 



}. 



w/zere Z is anmxc matrix, tr zs £fte £race on £fte boundary 
dM and v is the unit normal on dM. The subbundle V is 
a Dirac structure with respect to the bilinear form ((, )). 

PROOF. Using the fact that (Zf r , e x ) L 2 ( M ) = 
(/ r ,Z T e:c ) L ^ (M) forany/ r e C~(M) ande, G^fM), 
and applying the integration by parts formula 

(grad e x ,e d ) L 2 m (M) =(e x , div e d ) L ^ (M) 

+ (tie x ,tTe d • ^L^aMjj 

for e Cff(M), e d e C™ ( M), s imilar to Theorem 1 in 
van der Schaft & Maschke ( 2QQ2[ ), it is routine to show 
that V = V ± . □ 



In order to obtain a port-Hamiltonian formulation of 
reaction-diffusion systems, we start from the Gibb's free 
energy associated to the reaction system given by 



G[x) = x T Ln (^) + (x* - x) T l r 



(5) 



where x* is an equilibrium of the reaction network and 
l m denotes a vector of dimension m with all ones. It 
can be immediately checked that the gradient of G is 
the chemical potential /i, i.e., f^(^) = Ln (^-) =: ji{x). 
The energy, the Hamiltonian, of the reaction-diffusion 
system is Q = f M Gd£. 

Starting from the Dirac structure V given in Theorem 3, 
define the energy storage relations 



fx 



dx 
dt : 



dG 

dx 



(x). 



(6) 



The power-dissipation of the reaction part is 

f r = -£/C(x*)£ T Exp(e r ) , (7) 

where the mapping BJC(x*)B T Ex.p(-) satisfies 

e T r f r < for all e r e £ r • (8) 

This is due to the fact that the exponential function is 
strictly increasing, so 

r 

J=l 

• (exp (j v . (x)) - exp (75. (a?))) 

> 

for Kj(x*) > 0, j = 1, . . . ,r. 

Ordinarily, diffusion is treated as power-dissipation by 
thermal motion of particles and as such quantitatively 
is characterized by the diffusion matrix D introduced in 
(1), which explicitly refers to the state variable x. In the 
work at hand diffusion is modeled by termination of the 
diffusion port as 



ed = -Kd(fd) , 



(9) 



where !Z d : T d — > £d is in general a nonlinear mapping 
satisfying dissipation inequality 

e T J d < for all f d ^T d . 

The operator 1Z dl instead of acting upon the gradient of 
the state x, for the argument takes the gradient of the 
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co-energy variables (i.e., the chemical potential /i). We 
call this operator the energy diffusion operator. 

When the diffusion operator 7Zd is a matrix func- 
tion of the state x, these constitutive relations define 
the reaction-diffusion system in the port-Hamiltonian 
framework as 



9x , . /_ ✓ \ , dG, v 
— =divl 7Cd(x)grad— (x) 

- Z5/C(x*)5 T Exp ( Z T ^-(x) 

\ dx j no) 

/ 6 = ^ d (x)grad -q^(x) ■ ^|aM , 
with the Hamiltonian is (5 = J M Gd£ and G given by (5) . 



Because gradLn(^) = diag . . . , ^— ^ gradx, 

the system (10) is in the form (1) with 7Zd(x) = 
diag (#1, . . . , x m ) -D(x) and the reaction dynamics 
f{x) = -ZBJC(x*)B T Exp (Z T Ln (^)). 

Standard Model. The dynamical analysis of the bal- 
anc ed rea ction networks presented in | van der Schaft e£| 
al. ( 2Q12[ ) is given on the state space modulo the space 
oiequilibrium points. For the sake of thermodynamical 
consistency, we rewrite the system (10) into the form 
given in terms of the disagreement vector as 



^=div(^(x)grad(^))+/(x) 

65 = Ln (^) |aM 

fb = Rd(x)grad(J^ • v\qm , 



(11) 



where R d (x) := n d (x)diag (j^, . . . , and / is given 
by the right-hand side of (2) . 

The existence of solutions for t he systems (10) and (11 



is a comp l ex iss ue. The papers [Morgan| ( [l991| ); [Fitzgib-| 
bon et al. ( 1997[ ) do provide a working framework for the 
systems with sepa rable Lyapunov functio ns. Further- 
more, according to |Fitzgibbon et ~al ( |1997 ), the system 
does not generate spatial patterns. The problem of the 
existence of classical solution and spatial uniformity of 
the steady state in the presence of Neumann's boundary 
conditions for the semilinear reaction-diffusion systems 
we plan to address in a separate contribution. 

Passivity. Define the complex affinity as j(x) = 
Z T Ln (^r). Assuming the existence of a classical solu- 
tion to (10), as an immediate consequence we obtain 



the following energy balance 

= - (Z^(x),BJC(x*)B T E X p(Z^(x))) LU M) 
+ {fi(x),div(Tl d (x)gra,dfi{x))) L2m(M) 

= - (7(x), J B/C(x*) J B T Exp( 7 (x))) i \ (M) 
- (gTadn(x),TZ d (x)gTadn(x)} L2m{M) 
+ (li(x),n d (x)gradn(x) -v) L ^ dM) . 

Because the exponential function is strictly increasing 
the following inequality holds 

7 T B/CB T Exp(7) > 

for Kj(x*) > 0, j = 1, . . . , r, which immediately implies 

{Z T ii(x),B)C(x*)B T Ex V (Z^(x))) LUM) > 0. 

Furthermore, since 

(grad/i(x),^ d (x)grad/i(x)) L ^ (M) > 0, 

the passivity property holds 
d 



dt 



G < (e&,/&) L 2 



(dM) 



(12) 



which means that the Hamiltonian functional Q is non- 
increasing along solution trajectories of (10). The equal 
conclusion, of course, also holds for the system (11). 



5 Structure-Preserving Discretization 

A Single Species System. Firstly, let us consider the 
single component reaction-diffusion system 



dx 

— = div (D(x) grad x) + g(x) 

eb = AdM 

f b gradx • v\ dM , 



(13) 



where x,g,D e C°°(M), D(x) > for all x, and e b e 
C°°(dM) and f b e C°°(dM). 

In the framework of exterior geometry, we can identify 
0- and n- forms with scalar valued functions, while 1- 
and (n — 1) -forms are identified with proxy fields. This 
allows for the operators grad and div to be rewritten in 
terms of the exterior derivative d and the Hodge star * 
formally as: grad = d and div = — * d*; for more details 
see, e.g.,|Arnold et al\ (12010). 
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Following the exposition of Seslija et al. (2012b), let K 
be a homological simplicial complex obtained by trian- 
gulation of the manifold M. Assuming that K is well- 
centered, its circumcentric dual is *K = x *\>K, 
where *{K is the inte rior dual and *\>K is th e boundary 
dual, as explained in Seslija et al. ( 2Q12b|a ). 



The discrete analogue of an oriented manifold is an ori- 
ented simplicial complex, while differential forms are dis- 
cretized as cochains. A /c-cochain is a real- valued func- 
tion on the fc-simplices of K, which we will also call a 
discrete /c-form. Analogously, we define the space of dis- 
crete forms on *[K and *x>K. By f^(iT), f^(*iiT), and 
QjfoK) we denote the space of the primal /c-cochains, 
the dual /c-cochains, and the boundary dual /c-cochains, 
respectively. 

The discrete analogue of (13) is 



(*o)" 1 (dr- 1 *i D d (x)d°x + d£" 1 f b ) + g(x) 



dx 
dt 

e b = tr° x, 



(14) 



where the state x now lives on the set of verices of 
K, that is, x G ft° d (K), the input f b G ^ _1 (* b ^0, 
and the output e b G Q d (K). The positive-definite (dis- 



crete) diffusion matrix is x \-> Dd(x) G 



p N e x N e 



with 



dimQj (iT), while the opera tors *p, *i, d°, dp 



dp -1 , and tr° have been defined in Seslija et al. (2012b). 



The operator d° : Q(K) -> ft 1 (If) is nothing but the 
transpose of the incidence matrix of the primal skeleton 
(from the primal edges to the primal vertices). Further- 
more, d° = - (d[ l - 1 ) T and d^ 1 = (tr°) T . The discrete 
Hodge operator *i : ^(K) — » r2 n_1 (*ii\") is a diagonal 
matrix with the k-th entry being equal | *i a\ \ / \ a\ | , where 
o\ is the primal edge with the dual *\o~\. The matrix *o 
is a diagonal matrix whose /c-th element is | *i &®\/\a®\ . 

Remark 1 In fact, the model (14) slightly, but cru- 
cially, differs from the standard compartmental model 
on graphs, w here the matrix (*o) _1 does not appear 
A rca k\ $2011]/ . This implies that in the standard graph 
model *o = In? N = dim ft^(K) , that is, | *i a®\ = 1 
for all k = 1, . . . , Nj meaning that all the compartments 
are of equal volume. This fact does not surprise, since 
the graph formulation does not capture the geometric 
content of the underlying model. 

Multicomponent System. Let us now consider the 
reaction-diffusion system with m components (cf. 10). 
To each node of the primal mesh we associate reaction 
dynamics. That is, to a node a® we associate the state 
x 3 £ The geometric dual of a®, * s the dual 

volume cell which represents the j-th compartment (see 
Figure 2). The number of the compartments is N = 
dim ^ (if) = dim (*{!{). The compartments interact 
with each other through the diffusion modeled as follows. 




Fig. 2. A simplicial complex K consists of two triangles. 
The dual edges introduced by the circumcentric subdivision 
are shown dotted. The state vector x j = (xi,...,x m ) T is 
associated to the vertex Vj for each j G {l,...,iV}. The 
number of compartments for this example is N = 4. The 
shaded region, the dual cell *ii>2 of the vertex V2, represents 
the compartment with the state x 2 . 

By X denote the concatenated vector 

X=((x 1 ) T ) ...,(^) T ) T , (15) 

where x^ G Rip, and let 

F(X)=(f(x')\...j(x N ) T y (16) 

be the vector field which describes the reaction dynamics 
of all compartments, with the reaction kinetics f(x^) = 

-Z5/C(x*)5 T Exp (z T Ln (fQ) , j = 1, . . . , N. 

The open compartmental model of the reaction- 
diffusion system (10) is given by 



X -- 

e b = (tr <g> I m ) 



((*o) 1 ®I m )\&d 
X 

X*' 



X 

X* 



(tr ® I m y f b )+F(X) 



where (g) represents the Kronecker product, I m is the 
identity matrix of dimension m x m, R d (X) > aI m N e1 
a > 0, and N e is the number of edges of the primal 
mesh. The Laplacian matrix of the simplicial complex 
is A d = (d <g> I m ) T (*i <g> I m ) R d (X) (d <g> I m ). Note that 
we have used d to denote d° : 



(dp 



-l \ T 



and jj^r 



(cc* ) ' (cc* ) '•••'( x * ) 



The total energy of the system, the sum of energies of 
all compartments, is 



N 



G d (X) = J2Gi(xi)V a o, 

where is the vertex corresponding to the state x J and 
V a o is the n-dimensional support volume obtained by 

3 

taking the convex hull of the simplex a® and and its dual 
cell*i<7?. Since V a o = |a°||* i( 7°| = \^a% j = l,...,N, 
the total energy can be written as 



G d (X)-- 



N 
3 = 1 



(G 1 ,...,G N )* 1 N . (17) 
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The distributed chemical potential as the gradient of 
(17) is given as 



dGd 
OX 



: ■ /, I " ( A I ) • (18) 



Compartmental Model. Imposing zero- flux boundary 
conditions, /b = 0, leads to the closed compartmental 
model 

X= F*{X) =- ((^o)" 1 ® 7 m ) A d ^ + F(X), (19) 



with a positive initial condition X(0) = Xq G 



A simple but crucial observation is that F£ (X) > when 
X k = for any k = 1, . . . , raX and X G RT^. We have 



F£(X) = - Z^/C(x*)5 T Exp Z T Ln 



X*' 



where k = ((j — l)m + z) and X& = xj, while Z z is the 
z-th row vector of Z and (d <8) 7 m )j[ is the fc-th column 
of(d(g)J m ). 

When Xk — x\ = 0, the terms corresponding to the 
positive z-th diagonal element of the weighted Laplacian 
matrix BJC(x*)B T are all zero, while there is at least one 
term corresponding to a non-zero, and therefore strictly 
negative, off-diagonal element of B)C(x*)B T . This im- 
plies that < 0. Similarly, the matrix (*i (g) I m ) R d {X) 
is a diagonal strictly positive definite matrix, thus the 
terms corresponding to the positive ((j — l)m + i)-th 
diagonal element of the augmented Laplacian matrix 
(d / m ) T (*i0 Im) Rd{X) (d I m ) are all zero, while 
there is at least one off-diagonal negative element. Thus, 
S% < 0, and therefore F*{X) > when X k = 0. 

Lemma 4 Suppose that X : [0, t*] R™ N is any solu- 
tion of (19). Then for any k = 1, . . . ,mN: 

X fc (0)>0 => X k (t*)>0. 

PROO F. The proof is a repetition of the arguments 
given in 



Sontag (2001), Lemma 7.1, for a different class 



of systems. 

Suppose that k is so that X k (0) > 0. Let $ : R 2 R be 
the function which for t G [0, t*] and y G R coincides with 
17) :=^ fc * . . j/, X fc+1 (t),. . .,X mN (t)) 

and has <£(£, y) = $(0, ?/) for t < and y) = $(t*,y) 
for t > t*. Since F*(X) > when X k = 0, $(t,0) > 
for all t. For £ G [0, the scalar function := X k (t) 



satisfies y(t) = &(t,y(t)). We need to prove that y never 
vanishes. To this end, let ^(t,p) := $(t,p) - $(t,0), 
and let = with z(0) = 2/(0). 

Because \I/ is locally Lipschitz and is an equilibrium of 
z = $(t,z),z(t) > for alH. Furthermore, z = SSf(t,z) < 
$(t,z) for all t, and thus by comparison z(t) < y(t). 
Since y(t*) is well-defined, z(t) remains bounded, and 
thus is defined for t = t* . Hence, y(t*) > z(t*)> 0. □ 

Corollary 5 For the system (19) the positive orthant 
M.™ N is forward invariant. 

Jus like in the case of the system (2), in order to exclude 
the existence of possible boundary equilibria, we shall 
assume the global persistency property. 

Conjecture 6 Given X G R™^, all the trajectories 
t H> X(t) of (19) satisfy: liminf^ooX^) > 0. 

In the absence of the diffusion terms, the dynamics of 
the spatially discrete systems are decoupled, and as such 
coincide with the dynamics of the balanced reaction sys- 
tem (2). In this scenario all the compartments exhibit 
asymptotically stable dynamics, but the steady states of 
all the compartments, in general, are not identical. The 
following theorem shows that the compartmental model 
(19) is asymptotically stable with the spatially uniform 
steady state. 

Theorem 7 Consider the compartmental model of bal- 
anced mass action reaction network given by (19). For 
every initial condition X (0) G R™ 7 ^ the species concen- 
trations x 1 , . . . , x N as t — » oo converge to x 1 = • • • = 
x N e£. 



PROOF. In |van der Schaft et al.\ ( [2012D the authors 
have shown that G in (5) satisfies G{x ) = and G(x) > 
0, Vx ^ x*, and for every real c > the set {x G R+ | 
G{x) < c} is compact. This easily can be checked. Let Xi 
and x* denote the i-th elements of x and x* respectively. 
From the strict concavity of the logarithmic function 
z — 1 > ln(z), \/z G R+, with equality if and only if 

= 1. Putting z = ^-,we have x\ — X{ + Xi In ( ) > 0, 

with equality if and only if xi = x*. This implies that 

G(x) = YT=i (x* -Xi+Xi In (jfi)^ > 0, with equality 

if and only if Xi = x|, i = 1, • • • , m. Thus G has a strict 
minimum at x — x* and and G(x) > 0, Vx ^ x*. 

The above stated properties of G immediately imply that 
the function X \-> Gd(X) in (17) satisfies 

G d (X*) = 0, G d (X)>0, VX^X*, (20) 

and is proper, i.e., for every real C > the set {X G 



iN 



Gd(X) < C} is compact. 



In what follows we will show that G d (X) = ^^(X)X 
^(X) satisfies 



G d (X) < for all X G R^ 



mN 



(21) 
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and 



G d (X) = if and only if x 1 



x N £ £. (22) 



We look for the time derivative of the total energy: 



G d = = ((* ® I m ) Ln | ] X 



-Ln 



-J (*0 ®-?"m) T ((*o) _1 ®im) 



(23) 



• (d 7 m ) T (*! J m ) ^(X) (d / m ) — 
Ln (^) (*o^^) T ^(X). 



Since (* i m ) T ((*o) 1 ® ^m) = ^mTV, we have 
-Ln [ (d 7 m ) T (*! <g> J m ) iJ d (X) (d I n 

Lll( ^) (*0^m)V(I) 
= - ( (d (g) / m ) Ln (^)) (*1 ® 7 m) R d (X) (d ^ / m ) ^ 



X 

X* 



X 



N 



+Ei*^i^r/(* i ) 

i=l s v ' 



(d ® J m ) Ln (£) , R d (X) (d ® / m ) ^ 



iV 



£|*i<X°M^) 



1=1 



We compute the expre ssion along the lines of van 

der Schaft etal] ( |2Q12[ ) ? as 



= -/i T (x i )^/C(x*)5 T Exp(Z T /i(x i )) 
= - 7 T (x i )5/C(x*)5 T Exp( 7 (x i )) 

= E^i(7 < s,(^)-7P i (^))^(^*) 

• (exp (j Pj (x*))- exp (7^. (a: 4 ))) 

<0, 



since Kj(x*) > for j = 1, . . . , r, and the exponential 
function is strictly increasing. The summand in the third 



line of (24) is zero only if 75. (x*) — 77?. (#*) = for 
every j. This is equivalent to having B T j(x) = 0. Thus, 

^(2*) = only if 5 T 7 = £ T Z T Ln = 0. It follows 

that 

e R {x l ) =0 if and only if x z e £ (25) 
for all z = l,...,iV. 

For the contribution of the compartmental diffusion dy- 
namics we have 



dG X 

(d Im) Qx' R d( X ) ( d ® Im) 



^ ( ^ n ( T * 

k=l 



Ln 



a\al\ 



x" x J 



>0, 



where 7V e is the number of edges of the primal mesh, 
and x % and x J are states associated to the nods i and j, 
and k is the edge between nods i and j. Because ln(-) is 

an increasing function, ^Ln — Ln possesses 
the same sign as — , and hence > 0. Further- 



ed = if and only if x 1 



(26) 



Now, G d = if and only if = and sd = 0. The 
intersection of the two conditions (26) and (25) gives 
(22). 

Since Gd is proper (in W£ N ) and the state trajectory 
t \-> X(t) remains in R™ N , (21) implies that t \-> X(t) is 
bounded in R™^. Therefore, boundedness of t \-> X(t), 
together with equations (21) and (22), by LaSalle's in- 
variance principle imply that all the species concentra- 
tions x 1 , . . . , x N converge to an element in £. □ 

Remark 2 Theorem 7 remains unaltered if we replace 
the reaction vector field (2) by any vector function of the 
form (not corresponding anymore to mass action kinet- 
ics) 



f{pP) = -ZBK z (x^x*W$ (^ T Ln ^— J j , (27) 
( 24 ) where JC g (x^x*) := diag(^f (x^ x*), • • • ,K,f(x^x*)) > 



for all x 3 e 



and $ 



is a mapping 



®(yir - ,Vc) = diag(/i (yi ),..., f c (y c )), with the func- 
tions i = 1, . . . , c, all monotonically increasing. 

In fact, a recent paper Rao et al. ( 2Q13| ) shows that, for 
instance, Michaelis-Menten kinetics are in the form (27), 
where K g (x, x*) is a rational but strictly positive definite 
matrix. 
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6 Chemical Example 

We illustrate our analysis on a simple chemical reaction 
model 



x 1 + x 2 



:X 3 



:Xi + X4, 



(28) 



where X\ is enzyme, X 2 substrate, X 3 intermediate 
product, and X4 product. The first (binding) and third 
(unbinding) steps are reversible. Many reactions in the 
glycolysis metabolic pathway are of this type. 

For instance, glucose-6-phosphate isomerase (alterna- 
tively known as phosphoglucose isomerase or phospho- 
hexose isomerase) is an enzyme that catalyzes the con- 
version of glucose 6-phosphate (G6P) into fructose 6- 
phosphate (F6P) in the second step of glycolysis. The 
change in structure is an isomerization, in which the G6P 
(Xi) has been converted to F6P (X4). The freely re- 
versible reaction requires an enzyme X2 , phosphohexose 
iso merase , to proceed; for more details see, e.g., |Berg e1\ 



al. 



(2007) 



The dynamical model of (28) governed by mass action 
kinetics is given by 



£1 = -k{ orw Xl x 2 + (/4 orw + k[ ev )x 3 - k T 2 eY x x x A 
x 2 = -k{° rw x lX2 + k[ ev x 3 



x 3 = kf™ Xl x 2 - (k\ 
X4 



+ /4 orw )x 3 + k r 2 ev Xl x A 



(29) 



kl orw x 3 



k 2 X\X^. . 



It can be easily checked that x\ = x\ = 1, x\ = 
k[ orw /k{ ev , x\ = k[ orw kl orw /(kl ev k r 2 ev ) is one of the equi- 
libria of the system (29). The complex stoichiometric 
matrix Z, the incidence matrix £?, and the stoichiometr- 
ric matrix S for the reaction network (28) are 



/l0l\ 
1 
1 

V 001 / 



B- 



(l o\ 

-1 1 

V° -v 



S = ZB 



/-I l\ 
-1 

1 -1 

1 



Since Z Sl = (1,1,0,0) T , Z Vl = Z$ 2 = (0,0,1,0) T , and 
Z-p 2 = (1,0,0, 1) T , for the chosen x*, the diagonal bal- 
anced reaction constants (cf. Section 2) are 



^,forw 











The system (29) now can be rewritten into the form (2), 
while the dynamics under the influence of diffusion is 
given by the reaction-diffusion model (10). 




Fig. 3. Solutions of (29), in the presence of the diag- 
onal diffusion term diag(di, d 2 , <fe, d^Ax, on the one-di- 
mensional spatial domain M — [0, 1] with initial con- 
ditions xi(f,0) = 4£ + 0.3, x 2 (f,0) = 1.3£ 2 + 0.1, 
x 3 (£, 0) = 2 sin^ (0 + 0.2^ + 0.2, x 4 (£,0) = 3f + 0.1, and Neu- 
mann's boundary conditions. Diffusion coefficients are set to 
be di = 0.33, d 2 = 0.72, d 3 = 0.91, and d 4 = 0.67. The reac- 
tion rates are ^ orw = 0.1, k{ ev = 0.4, /4° rw = 0.3, k r 2 ev = 0.5. 
Upon the transient phase the system reaches a steady state 
x** = x(f,oo) = (2.1856, 1.7557, 0.9602, 0.2638) T uniform in 
space. Immediately, we verify that x** is a thermodynamical 
equilibrium, S T Ln(x**) = S T Ln(x*). The number of com- 
partments used in N = 20. 

The spatially uniform asymptotic behavior predicted by 
Theorem 7 is demonstrated with the simulation in Fig- 
ure 3. The elements of the matrix in the system 
(19) given in terms of the standard diffusion matrix D 
are Rd := (diag(^r, . . . , Ar)D) (g) I m , where in our case 

m — 4. Eliminating the effects of diffusion leads to a 
spatially nonuniform steady state. 



7 Concluding Observations 

We have provided a geometric formulation of reaction- 
diffusion systems with a thermodynamical equilibrium. 
Diffusion as an isolated process is associated with a 
homogenizing effect that eliminates the gradients of the 
constituents and eventually leads to uniform spatial 
state. However, diffusion in combination with reaction 
dynamics can produce spatially heterogenous patterns. 
We envision that the Hamiltonian perspective pre- 
sented in this paper will utilize the systems analysis 
of reaction-diffusion system and foster new insights in 
compartment al systems design. 

Control for reaction-diffusion systems in the port- 
Hamiltonian framework can be understood as the cou- 
pling of a reaction-diffusion system to an additional 
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port-Hamiltonian system that plays the role of the 
controller. This, among others, enables the applica- 
tion of passivity base techniques in control synthesis 
for react ion- diffusion systems. Only by having accurate 
structured discretization can one hope to approach this 
challenging enterprise. 

A model obtained by structure-preserving scheme for the 
spatial discretization is a compartmental model, which 
exhibits a striking similarity with con sensus dynamics 
|01fati-Saber et a/.| ( |2QQ7| ); |Cortes| ( |2QQ8[ ). Exploring this 
resemblance is a very appealing research direction. 
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